Intermittency in Two-Dimensional Turbulence with Drag 
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We consider the enstrophy cascade in forced two-dimensional turbulence with a linear drag force. 
In the presence of linear drag, the energy wavenumber spectrum drops with a power law faster than 
in the case without drag, and the vorticity field becomes intermittent, as shown by the anomalous 
scaling of the vorticity structure functions. Using previous theory, we compare numerical simula- 
tion results with predictions for the power law exponent of the energy wavenumber spectrum and 
the scaling exponents of the vorticity structure functions (,2 q - We also study, both by numerical 
experiment and theoretical analysis, the multifractal structure of the viscous enstrophy dissipation 
in terms of its Renyi dimension spectrum D q and singularity spectrum f(a). We derive a relation 
between D q and <^2 q , and discuss its relevance to a version of the refined similarity hypothesis. In 
addition, we obtain and compare theoretically and numerically derived results for the dependence 
on separation r of the probability distribution of Spuj, the difference between the vorticity at two 
points separated by a distance r. Our numerical simulations are done on a 4096 x 4096 grid. 



I. INTRODUCTION 



Two-dimensional Navier-Stokes turbulence has at- 
tracted much interest because of its relevance to a variety 
of natural flow phenomena. Examples are plasma in the 
equatorial ionosphere and the large-scale dynamics 
of the Earth's atmosphere and oceans Q. In the lab- 
oratory, experiments that are close to two-dimensional, 
such as soap film flow 0,0 and magnetically forced strat- 
ified flow [5|, have been conducted. In addition, rotating 
fluid systems [j| are used to study quasi-two-dimensional 
turbulence and its relevance to large-scale planetary flow. 
Ref. gives a review of some recent experiments in two- 
dimensional turbulence. 

In many of the situations involving two-dimensional 
turbulence, there are regimes where drag is an important 
physical effect. In the ionospheric case, there is drag 
friction of the plasma as it moves relative to the neutral 
gas background (due to ion-neutral collision). For geo- 
physical flows and rotating fluid experiments, viscosity 
and the no-slip boundary condition give rise to Ekman 
friction. In this case the three-dimensional flow is of- 
ten modeled as two-dimensional outside the layer adjoin- 
ing the no-slip boundary, and the effect of the boundary 
layer manifests itself as drag in the two-dimensional de- 
scription. For soap film and magnetically forced flows, 
friction is exerted on the fluids by surrounding gas and 
the bottom of the container, respectively. In all these 
cases, the drag force can be modeled as proportional to 
the two-dimensional fluid velocity v, thus the describing 
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Navier-Stokes momentum equation becomes, 
d'v 1 

-gj; + (v- V)u= --\7p + isV 2 v~ iiv + f , (1) 

where p is the fluid density, /i is the drag coefficient, v 
is the kinematic viscosity and / is an external forcing 
term. In two-dimensions, the system can be described 
by the scalar vorticity field u> whose equation of motion 
is obtained by taking the curl of Eq. QJ, 



duj _ 2 

— h V ■ Vld = V\ U) — IIU! + f u 

at 



(2) 



with uj = z ■ V x v and f u — z ■ V x /, z being the unit vec- 
tor perpendicular to the plane. In our studies, the forcing 
will be taken to be localized at small wavenumbers (k) 
with a characteristic wavenumber kf, and incompressibil- 
ity, V • v = 0, will be assumed. 

According to Kraichnan, for two-dimensional turbu- 
lence with no drag and very small viscosity, for k 3> fc/, 
up to the viscous cutoff kd, enstrophy cascades from 
small to large k ,8j. As a result, the energy wavenum- 
ber spectrum E(k) has a power law behavior with log- 
arithmic correction 0, E{k) ~ fc- 3 [ln(fc/fc / )]- 1 /3 for 
kf <C k kd- In the presence of a linear drag, the 
energy spectrum drops with a power law faster than the 
case with no drag [l(J, EH [12 > 



E(k) ~ fc-( 3+ «> (£ > 0) 



(3) 
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and there is no logarithmic correction. A result similar 
to Eq. ||2J has been obtained for the closely related prob- 
lem of chaotically advected finite lifetime passive scalars 
[HHlHil. (See SectionHHlfor a discussion of the rela- 
tionship between vorticity in two-dimensional turbulence 
with drag and finite lifetime passive scalars in Lagrangian 
chaotic flows.) 
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Furthermore, the vorticity field is intermittent as indi- 
cated by (i) anomalous scaling of the vorticity structure 
function, (ii) scale dependence of the vorticity difference 
distribution function, and (iii) multifractal properties of 
the enstrophy dissipation field. Intermittency in the 
closely related problem of chaotically advected finite life- 
time passive scalars was originally studied by Chertkov 
|13| using a model flow in which the velocity field was spa- 
tially linear and 5-correlated in time (white noise). This 
model successfully captures the generic intermittent as- 
pect of the problem [l(|. With respect to the situation 
of interest to us in the present paper, it is a priori un- 
clear to what extent the (5-correlated model yields results 
that approximate quantitative results for intermittency 
measures obtained from experiments (numerical or real) . 
This will be discussed further in Sees. IIII Bl and llVl 

The structure function of order 2q, S2 q (r) is defined 
as the (2q) th moment of the absolute value of the vor- 
ticity difference Spcu = u>(x + r) — u(x). Assuming that 
the system is homogeneous and isotropic, the structure 
functions depend onr= |r| only. For the case with drag 
(fx > 0), it is found that, in the enstrophy cascade range 
i^d 1 ^ r ^ kj 1 ), S2q(r) scales with r as, 



S 2q (r) = (l^n ~ r 



r (2q 



(4) 



with (^2q > 0. Furthermore, the vorticity structure func- 
tions show anomalous scaling; that is, C^q is a nonlinear 
function of q. The nonlinear dependence of C,2q on q in- 
dicates that the vorticity field is intermittent. In con- 
trast, in the absence of drag (/x = 0), it is predicted that 
w{x) wiggles rapidly {i.e., on the scale A;^ 1 ) and homo- 
geneously in space. In terms of Eq. (0J, this corresponds 

toc 2 g=ofor /x=o nmm. 

The intermittency of the vorticity field also manifests 
itself as a change in shape or form of the probability dis- 
tribution function of the vorticity difference 5?uj with the 
separating distance r. It can be shown that if the prob- 
ability distribution function P r (Xp) of the standardized 
vorticity difference, 



X^ 



(5) 



is independent of r, then & q increases linearly with q: 
From Eq. |g} and Eq. JSJ, 



S 2q (r) = (\SM 



) q J \X?\ 2q P r (X?) dXp 



(6) 



and, if P r (Xp) is independent of r, for r in the cas- 
cade range, then the only r-dependence of S2 q (r) comes 
through the term (\Spu>\ 2 ) q , implying C^ q = 9C2; that is, 
C,2 q is linearly proportional to q. Such collapse of P r (Xp) 
for different values of r has been observed in an experi- 
ment where drag is believed to be unimportant. When 
the effect of drag is not negligible, P r (X?) changes shape 
and develops exponential or stretched-exponential tails 
as r decreases, so that Eq. (jHJ admits nonlinear depen- 
dence of £29 on q. 



The intermittency of the vorticity field is closely re- 
lated to the multifractal structure of the viscous enstro- 
phy dissipation field. We define the enstrophy as iv 2 /2. 
From Eq. © , the time evolution of the enstrophy is then 
given by 



d 1 ur 
di \~2 



^ =r = -V- 



- u\Vw\ 



-2a*( y ) +UJ f^ 



(7) 



We identify the second term on the right hand side of 
Eq. as the local rate of viscous enstrophy dissipation 

V, 



(8) 



The multifractality of the viscous enstrophy dissipation 
can be quantified by the Renyi dimension spectrum of 
a measure based on the vorticity gradient. Imagine we 
divide the region TZ occupied by the fluid into a grid of 
square boxes of size e, we define the measure pi of the i th 
box IZi (e) as 



(9) 



The Renyi dimension spectrum |l9j based on this mea- 
sure is then given by, 



D„ 



-1- lim lim l ° g ^ 
— 1 e^o i/->o loge 



(10) 



The definition Eq. (|10|l was introduced in the context 
of natural measures occurring in dynamical systems by 
Grassberger [2(j, and Hentschel and Procaccia [2l|]. In 
the case with drag, we find that the dimension spectrum 
for pi is multifractal; that is, D q varies with g, in contrast 
to the case with no drag in which D q — 2, indicating that 
the measure is uniformly rough. 

The measure pi can also be described in terms of the 
singularity spectrum /(a) |22j| . In particular, to each box 
IZi, we associate a singularity index aj via 



a, 



logP» 
loge 



(11) 



and the number of boxes N{a)da with singularity index 
between a and a + da is then assumed to scale as 



N(a) 



(12) 



/(a) can loosely be interpreted as the dimension of the 
set of boxes with singularity index a [23|. When f(a) 
and D q are smooth functions, f(a) is related to D q by a 
Legendre transformation (2^ • The multifractal nature of 
the viscous enstrophy dissipation in the presence of drag 
implies that the /(a) spectrum, defined by Eq. (|llfl and 
Eq. I|12|l . is a nontrivial function of a. 
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The subject of this paper is the relation of P r (Xp), the 
exponents, £ and and the fractal dimension D q to 
the chaotic properties of the turbulent flow. In chaotic 
flows, the infinitesimal separation between two fluid par- 
ticle trajectories, Sx(t) typically diverges exponentially. 
The net rate of exponentiation over a time interval from 
to t for a trajectory starting at xq is given by the finite- 
time Lyapunov exponent, h defined as, 

, , , 1 , \Sx(t)\ 

h ^=i l ^\sm\- (13) 

At a particular time t, h(t;xo) in general depends on 
the initial positions xo and the initial orientation of the 
perturbation 5x(0). However, for large t, the results for 
h(t;xo) is insensitive to the orientation of 6x(0) for typ- 
ical choices of 5x(0), and we neglect the dependence on 
Sx(0) in what follows. The distribution in the values of 
h for randomly chosen xq can be characterized by the 
conditional probability density function P(h \ t). In sub- 
sequent sections, we shall briefly review the theories that 
relate £ [Hj and C,2q jUJ to the distribution P(h \ t) and 
the drag coefficient /i. We then derive an expression for 
P r (Xp) and a relation between D q and C,2q- We apply 
these theories to compute £, (,2 q , Pr(X?) and D q for tur- 
bulent flows governed by Eq. (|2J) and the theoretical re- 
sults are compared to those obtained from direct numer- 
ical simulations. 

We perform our simulations on a square domain of size 
[—it, 7r] x [— 7T, 7r] with periodic boundary conditions in 
both directions. The viscous term in Eq. @ is replaced 
by a hyperviscous damping — j/V 8 w with v — 7.5 x 1CP 25 
and fu>(x,y) = cos2x is used for the source function 
of the vorticity. Our use of a hyperviscosity is simi- 
lar to what is often done in numerical studies of three- 
dimensional turbulence and, for a given numerical reso- 
lution, has the desirable effect of increasing the scaling 
range where dissipation can be neglected, while, at the 
same time it is hoped that this change in the dissipa- 
tion does not influence the scaling range physics. For 
wavenumber k < 6, we have fi = 0.1 and this provides 
an energy sink at the large scales. For k > 6, we will 
consider the cases of /1 = 0.1 and /1 = 0.2. As we shall 
see in Section Hi Al when drag is present, the large k vor- 
ticity components can be considered as being passively 
advected by the small-A; flow components. Applying the 
same small-A; drag (i.e., for k < 6) allows us to 

compare the effects of drag at small scales while keeping 
similar large scale dynamics of the flows. For all the nu- 
merical results presented here, we use a spatial grid of 
4096 x 4096 and a time step of 0.00025. Starting from 
random initial conditions for the vorticity field, Eq. @ 
is integrated u sing a time split-step method described in 
detail in Ref. |25|. The system appears to reach a sta- 
tistical steady-state after about 40 time units. FIG. Q] 
shows snapshots of the vorticity field for the cases of 
H(k > 6) = 0.1 and fx(k > 6) = 0.2. We note that the 
vorticity field for the case with a larger drag shows less 
fine structure and the large scale vortices tend to have 




FIG. 1: Snapshots of the vorticity field uj(x,y) at t = 61 
for the case fi(k > 6) = 0.1 (upper) and at t = 65 for the 
case fi(k > 6) = 0.2 (lower) (/x = 0.1 for k < 6 for both 
cases). Light areas are regions of large positive values of the 
vorticity, and dark areas are regions of negative vorticity of 
large magnitude. The scales used in the two plots are not the 
same. 



shorter lifetime. At any given moment, there are typi- 
cally 3 - 5 large vortex structures visible in the system. 
In the steady-state, vortices are continuously created and 
destroyed. 
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II. ROLE OF THE FINITE-TIME LYAPUNOV 
EXPONENT 

A. Passive Nature of the Small-Scale Vorticity 

The theory that we will use is based on the approxi- 
mation that the high k components of the vorticity field 
are passively advected by the large scale structures of 
the flow. This can be justified by the following argument 
given in Ref. ^3] • The Lyapunov exponent h of the flow is 
the mean rate of exponentiation of differential displace- 
ment Sx following the flow, where d(6x)/dt = Sx ■ Vv. 
Thus one can crudely estimate the Lyapunov exponent 
as 



h ~ (||Vv|| 2 )2 ~ ^jf k 2 E(k)dk , (14) 

where ||Vv|| 2 = (dv x /dx) 2 + (dv y /dy) 2 + {dv x /dy) 2 + 
(dv y / dx) 2 . Assuming the limit of infinite Reynolds num- 
ber and power law behavior of E(k) valid as k — ► oo, the 
integral in Eq. I|14|l diverges at the upper limit unless £ in 
Eq. © is positive. That is, the velocity field v is not dif- 
ferentiate (h and Vv are undefined) unless £ > (alter- 
natively, if £ < and viscosity imposes a cutoff to power 
law behavior of E(k) at k ~ kd, then, although ||Vu|| 2 
is now finite, the integral in Eq. (|14|l is dominated by 
velocity components at the shortest wavelength). From 
Eq. ifHjl. for £ > 0, we have h ~ kj ^ 2 . This means that 
h, which characterizes small scale stretching, is deter- 
mined by the largest scale flow components. Since £ > 
in the case where drag is present, Vv is predominantly 
determined by the largest spatial scales. Thus the Lya- 
punov exponents provide information on the evolution 
of the distance between fluid elements whose separation 
is finite but small compared to kj 1 . This will allow us 
to approximate the evolution of vorticity field compo- 
nents with wavenumbers in the range kf <C k < kd us- 
ing Lyapunov exponents that result primarily from the 
large spatial scales k <~ kt. That is, the vorticity field 
at wavenumber kf <C k < kd evolves via approximately 
passive advection by the large scale flow. (Note that for 
£ < such an approach would not be applicable since the 
Lyapunov numbers would provide an estimate of separa- 
tion evolution only for distances less that kj 1 which is 
outside the dissipationless power law range.) 

Note that the case without drag corresponds to £ — 0, 
which is marginal in the sense that it is on the border- 
line of the condition for differentiability of the velocity 
field. In other situations of marginality (e.g., in critical 
phenomena), it is often found that there are logarithmic 
corrections to power-law scaling, and this may be thought 
of as the origin of Kraichnan's logarithmic correction to 
the A; -3 enstrophy cascade scaling of E(k). 
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FIG. 2: Probability distribution function of the finite-time 
Lyapunov exponent P(h \ t) for different time t for the case 
p = 0.2. 



B. Distribution of Finite-time Lyapunov Exponent 

As mentioned earlier, the exponential divergence of 
nearby trajectories in a chaotic flow over a time interval 
to t can be quantified by a finite-time Lyapunov exponent 
h(t;xo) defined in Eq. JT3J. In the limit t — > oo, h(t;xo) 
will approach the usual infinite-time Lyapunov exponent 
h for almost all initial conditions Xq and almost all initial 
orientations of Sx. At finite time, the dependence of h 
on x results in a distribution in the values of h which 
can be characterized by the probability density function 
P(h 1 1). That is, if xq is chosen randomly with uniform 
density in the region of the fluid flow, and if the initial 
orientation of Sx is arbitrarily chosen, then we can define 
a probability distribution P{h \ t) such that P(h | t)dh is 
the probability that h(t; xo) lies between h and h+dh. As 
t increases, P{h \ t) will become more and more sharply 
peaked at h and approach a delta- function at h as t — > oo. 

The passive nature of the small-scale flow allows us to 
use the procedures described in Ref. [2£j to obtain his- 
togram approximations to P(h 1 1). FIG.[5]shows P(h 1 1) 
at different t for the case of /z = 0.2. As time increases, 
P(h 1 1) becomes sharply peak around a positive value 
of h showing that the flow is chaotic. In particular, for 
t = 20, P(h 1 1) has its peak at h « 0.20. The function 
P(h 1 1) shows similar behavior for the case fi = 0.1 with 
a peak occurring at h « 0.24 for t = 20. The smaller 
damping in the latter case apparently allows the fluid to 
undergo bigger stretching in a given amount of time. 

Based on the argument that h(t; xq) can be considered 
as an average over many independent random realiza- 
tions, P(h 1 1) is approximated by the following asymp- 
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h 

FIG. 3: G(h) (see Eq. 1151 1 for different time t for the case 
jj,(k > 6) = 0.2. The inset shows G(h) obtained at t = 20 
(solid line) and the quadratic approximation to G(h) obtained 
by the method described at the end of Sec. 1111 fjl fdas tied line). 



totic form ( 23] and references therein), 

P{h\t)~ } J^j£e-« s W (15) 

for large t. Eq. (|15|l has been shown to be true for the 
generalized baker's map |2("il| and numerically verified for 
cases where there are no KAM surfaces, for example, see 
Ref. . The function G(h) is concave upward, G"{h) > 
0. It has the minimum value zero, occurring at h = h, 
G'(h) = G(h) = 0. The significance of Eq. JTjj is that it 
expresses a function of two variables P(h \ t) in terms of a 
function of a single variable G(h). Eq. (|15|l will be used 
in the development of the theory. We note that G(h) 
is completely specified by the flow v(x, t) and hence, is 
dependent on the value of fi. 

The function G(h) can be approximated at large t from 
P(h 1 1) using the following relation, 

G(h) ss K- iln P(h 1 1) , (16) 

where K is determined by the condition that the min- 
imum of G{h) equals zero. FIG. shows the G(h) ob- 
tained from the corresponding P(h \ t) shown in FIG. 
The large number (4096 2 ) of fluid trajectories used in 
the generation of each P(h 1 1) allows us to obtain G{h) 
for a large range of h. As can be seen in FIG. |3 for 
large enough t, the graphs of G(h) for different t more or 
less collapse onto each other, showing that Eq. (|15|) is a 
good approximation to P(h \ t) for the flows we consid- 
ered. G{h) shows similar behavior in the case /z = 0.1. 



III. THE EXPONENTS £ AND £2, 
A. Review of Theory 

We consider the scaling of the vorticity structure func- 
tions in the limit v — > + . The results have previously 
been given in Ref. Hp. 123 , which treats finite-lifetime 
passive scalars 0, Il4 Il5| rather than vorticity in two- 
dimensional turbulence with drag, and in Ref. |28| for 
two-dimensional turbulence with drag. Here, we re- 
derive the previous results (our derivation is different 
from that in Ref. [HI . l2~i| and more complete than that 
in Ref. |H). 

From Eq. f2j) with initial condition ui(x, 0) = 0, we 
have Hi, 

Jo 

Jt-r(r) 

+ f SUxit'^'-^dt' . (17) 
Jo 

where = fu[x (f) +5x{t')] - f u [x(t% x(t') and 

x(t') + 5x(t') are the two trajectories that pass through 
x and x + f at t' = t; thus 5x(t) = r. Since the points 
x and x + r are specified at time t, chaos in the flow 
implies that the backward evolved trajectories to times 
t' < t diverge exponentially from each other, Sx(t') ~ 
fe h ^ t \ until |<5x(t')| reaches the system size L, past 
which \5x(t')\ remains of order L. We define r(r) such 
that \Sx(t — T)\ ~ L ~ re hT . We have linearized Sf u [x(t')] 
for t — r(r) < t' < t when the two trajectories are close 
together. On the other hand, for < t' < t — r(r), 
\8x(t')\ ~ L and 5f,jj\x(t')\ thus fluctuates with roughly 
constant amplitude, 5f u ~ f u . Thus the integral from 
to t - r(r) in Eq. (HJ is of the order of e~^ T p^ . 
The integral from £ — r(r) to t in Eq. 117(1 shows the 
same ~ e~^ T behavior if h > \x but is of the order of 
e~ hT if h < fx. Hence, using the definition of S2 Q {r) and 
replacing the average over x by an average over r at fixed 
exponentiation A = ln(L/r) = /it, we find 

/>— poo 

S 2q (r)~ / dr R(t I A)e- 29 ' IT + / dr R(r | A) e - 2?A . 

(18) 

The conditional probability density R{t | A) of r for a 
fixed A is related to P(h \ t) by psj 

i?( T | A) « — / dhP(h\r) . (19) 
dr 

Using the asymptotic form Eq. l|15f) for P(/i | t), the in- 
tegral Eq. 1(18(1 is performed using the steepest descent 
method with A as the large parameter. Thus we obtain 
that the structure function scaling exponents, £29 defined 
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in Eq. Q, is given by 

&q = 

where 



mm{2q,H q (h)} , 

h 



H q (h) = 



G(h) + 2q(i 
h 



(20) 



(21) 



Equations ij^Ufran d (f^T|) yield the previous passive scalar 
result of Ref. [13| if G(h) is assumed parabolic, G(h) = 
(constant)(h — ft,) 2 , which is a consequence of the tem- 
porally white noise velocity field model of Ref. . 

We now consider the energy density e(k) which is given 
by 



e(k) 



1 \v x (k)\ 2 + \v y (k)\ 



(2nf 



2L 2 



(22) 



where v x {k) and v y (k) are Fourier transforms of the x and 
y components of the velocity field v(x, y). The wavenum- 
ber energy spectrum E(k) is then defined as 



E{k) = / dk'5(\k'\ - k)e{k') . 



(23) 



A previous theory pH . lllj | relates the energy wavenumber 
spectrum exponent, £ defined in Eq. © to G(h). The 
result is 



£ = mhi{#i(ft)} 

h 



(24) 



Thus, £ and ( 2q are related to the properties of the flow, 
namely the drag coefficient \i and the distribution of the 
finite-time Lyapunov exponent ft. 

Let ft = ft* be the value of ft at which H q {h) is min- 
imum. We now show that ft* increases with q. Letting 
(3 > a, by the definition of ft*, H' a {h* a ) = H' (h*p) = 0, 
it follows that G'(h%) - H p {h* ) = G'{h* a ) - H a (h* a ). 
Since by Eq. lf2TJl. Hp(h) > H a {h) for all ft, we have 
G'{h*p) > G'(h* a ) which implies h% > ft* due to the fact 
that G"(h) > for all ft. Moreover, putting a = gives 
ft* > ft for all q. 



B. Comparison of Theory and Numerical Results 

To apply the theory Eq. l(2T7)l and Eq. (|2*4*jl to numeri- 
cally determine the exponents £ 2q and £ of our turbulent 
flow governed by Eq. (J2J, we first let 



Q 2q = mm{H q (h)} . 

h 

Using Eq. (|2J , we re- write Eq. as 

m r\ — h — 



(25) 



(26) 



Since the function minimized in Eq. l|26[l has minimum 
value zero, we can multiply it by any positive function of 




6 8 10 12 14 16 18 20 

t 

FIG. 4: Log-linear plot of the partition function T{C$ q ,t) for 
fj, = 0.1. The dotted lines are linear fits. 



ft, and the minimum will still be zero. Thus for ft, > we 
can multiply the minimized function by ft to obtain 



min{G(ft) - h( 2q } = ~2q\i 



(27) 



Using Eq. l(To|) and Eq. (|2*7J) . we see that steepest descent 
evaluation of the following integral for large t yields 



P(h | t)e h ^" t dh 



(28) 



Thus we define the partition function |29| | 

T(z, t)= f e zh ^ 3a ^ dx , (29) 



(30) 



in terms of which Eq. l|28fl becomes 



r(Ca„*)~e 



2qfit 



We numerically compute the partition function Eq. I|29|l 
using the approximation, 



M 



M ^ 

" =i 



(31) 



employing many (Af=4096 2 ) spatially uniformly dis- 
tributed initial conditions xoi (i = 1, 2, . . . , M). 

For different values of £ 2 q> we then plot InTfaq, t) ver- 
sus t. FIG. ^ shows samples of \nT(( 2q ,t) for the case 
[i = 0.1. As expected, for large t, lnT(( 2q ,t) is linear 
in t and the slope, which can be estimated using a lin- 
ear fit, gives the corresponding value of q for each ( 2q . 
FIG.Elplots ( 2q versus q obtained in this way for fj, = 0.1 
and fi(k > 6) = 0.2 (open circles and squares) together 
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FIG. 5: as a function of q from the method of FIG. [I] 
for /i = 0.1 (circle) and > 6) = 0.2 (square). The solid 
lines are sixth degree polynomial fits to circles and squares. 
Results using Eq. (5) of Ref. |l3| are shown as dashed lines. 



what surprisingly, however, this does not lead to much 
difference in the numerical values of the structure func- 
tion exponents. This is shown in FIG.JSJwhere the results 
using Eq. (5) of Ref. are plotted as dashed curves. 
In fact, the difference is within the error of our numeri- 
cal experiments that directly determine the values of £29 ■ 
Thus for this case the parabolic approximation provides 
an adequate fit to the data, although this could not be 
definitely predicted a priori. 



1. Energy Wavenumber Spectrum 

As illustrated in FIG. 03 theoretical predictions for £, 
denoted as £th, are obtained using the method described 
above. The results are given in TABLE [I] To verify 
the theoretical results, we compute the energy spectrum 
directly from the numerical solution of Eq. J3J on a 4096 x 
4096 grid using Eq. l|23fl . which can be written in terms 
of the vorticity as 



E(k,t) = 




with the corresponding sixth degree polynomial fits (solid 
lines). By Eq. and Eq. (|2U|) . the exponents £ and 
(,2 q can then be determined from FIG. 03 by £ = (2 an d 
C-2q = min{2<7,C 29 }. 

It is also possible to compute £ and £2? directly from 
Eq. (|24() and Eq. (|20() using a fourth degree polynomial 
fit of the numerically obtained G(h). We find that the 
two methods give similar results, but that the method of 
Eq. lp?T|) is easier to compute. 

We also test the applicability of the result of Ref. 0] , 
which employs a model in which the velocity field is 5- 
correlated in time. This model may be shown to cor- 
respond to Eq. lt2"Ull and Eq. (|2T)> with G(h) parabolic, 
G(h) = a(h — h) 2 where h is the infinite time Lyapunov 
exponent. With this form of G[h), an explicit analytic 
expression for the structure function exponents can be 
obtained in terms of the constant o and h; see Eq. (5) 
of Ref. 0. In order to apply this result to a specific 
flow, we need to formulate a procedure for obtaining a 
reasonable value of a from the flow (h is well-defined and 
numerically accessible by standard technique). For this 
purpose, we note that, within the (5-correlated approxi- 
mation, the total exponentiation h(t, xo)t experienced by 
an infinitesimal vector originating at position xo under- 
goes a random walk with diffusion constant D = a/4. 
Thus we obtain a as one half of the long time slope of a 
plot of {(h - h) 2 )t 2 versus t. The inset to FIG. shows 
as a solid curve G(h) obtained using Eq. l(T(j|) (as previ- 
ously described) along with the parabolic model (dashed 
curve) with a determined by the above procedure. There 
appears to be a substantial difference in the important 
range h > h where the saddle points are located. Some- 



where Cj(k' \t) is the Fourier transform of tu(x,t). The 
time averaged energy spectrum E(k) is obtained by av- 
eraging E(k,t) at every 0.1 time unit from t = 41 to 
t = 75. FIG.Elshows a log-log plot of E(k) versus k for 
the two different values of /i(fc > 6) we considered. In 
both cases, a clear scaling range of more than a decade 
can be observed. We measure the scaling exponents by 
linearly fitting E(k) in the scaling range. The results, 
denoted £,dns, are shown in TABLE [I] Good agreement 
is found between the numerical and the theoretical re- 
sults. These results are consistent with those of previous 
work in and [28| which use grids of 1024 x 1024 and 
2048 x 2048, respectively. 



2. Vorticity Structure Functions 

Numerical tests of theoretical predictions of structure 
function exponents of finite-lifetime passive scalar fields 
advected by simple chaotic flows have been performed 
in Ref. [24| and |28j|. To test the analogous theoretical 
predictions, Eq. (|2U|) and Eq. (|2~T)l . in the case of two- 
dimensional turbulence with drag, we define the averaged 



TABLE I: Comparison of the values of 5 obtained from nu- 
merical simulations to the theoretical results. 



H(k > 6) 


C*h(— C2,i/i) 


£,DNS 


C2, BiVS 


C2,B(fc) 


0.1 


0.63 


0.61 


0.66 


0.68 


0.2 


1.10 


1.12 


1.16 


1.14 
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FIG. 6: Energy wavenumber spectra. The dotted lines are 
the corresponding linear fit in the scaling range. 



structure functions of order 2q as 



dr' 



Si q {r)= I ^- r 5{r-\r'\)(\8M 2q ) 



(33) 



The angled brackets denote average over the entire re- 
gion occupied by the fluid. The angular dependence of 
{\6f*aj\ 2q ) is averaged out in Eq. (|33|l by the integration 
over r'. Using Eq. (|33(l with oj(x,t) obtained from the 
numerical integration of Eq. @, we compute S2 q (r) from 
t = 41 to t = 75 at every 1 time unit and take the average 
of the results obtained. 

Following the scheme described above, we calculate the 
time-averaged structure functions for q ranging from 0.0 
to 2.0. FIGS.EJa) andlHfa) show samples of the results 
for the cases \x =0.1 and 0.2, respectively. The distance 
r is measured in the unit of grid size. For all values of 
q we have studied, the structure functions show a clear 
scaling range that is long enough to allow an estimate 
of the scaling exponents, &q- The scaling range of the 
structure functions in real space roughly corresponds to 
that of the energy spectrum in fc-space. The values of 
&q are obtained by measuring the slope of the structure 
functions in the scaling range using a linear fit. Results 
for (,2q/C,2 are shown as circles in FIGS. |7fb) andHJb). 
The measured values of (2, denoted as C2, £>ivs, are given 
in TABLE [I] We then obtain theoretical predictions to 
(21J from the polynomial fit of £25 shown in FIG. 03 follow- 
ing procedures described at the beginning of this section. 
The results are shown as crosses in FIG.[7{b) andOHb) for 
the two cases of /z we studied. The predicted values of £2, 
denoted as C2,th, are given in TABLE IJ The numerical 
and theoretical results agree reasonably well for all q's. 
In reference to Eq. (|20|l we note that, for all the cases 
we numerically tested, we found that Qi q — Qiq < 2q for 




q 

FIG. 7: For the case of /i = 0.1: (a) structure functions of 
vorticity difference, for various orders q between 0.1 and 2.0; 
the dotted lines are linear fits in the scaling range, (b) Plot of 
C,2 q /C,2 obtained from numerical simulations (circle) and from 
Eq. 1201 (cross), for different values of q; the solid lines are 
polynomial fits to the circles and the crosses (cf. Eq. (154H ). 



all q. This is consistent with the general result that for 
fi < h, C,2 q < H q (h) < 2q. We note, however, that this 
need not hold in general, particularly for large /1. 
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FIG. 8: For the case of /i = 0.2: (a) structure functions of 
vorticity difference, for various orders q between 0.1 and 2.0; 
the dotted lines are linear fits in the scaling range, (b) Plot of 
C29/C2 obtained from numerical simulations (circle) and from 
Eq. 1201 (cross), for different values of q; the solid lines are 
polynomial fits to the circles and the crosses (cf. Eq. 1541 1. 



C. Discussion 

1. Energy Wavenumber Spectrum 

The presence of drag gives a power law for the en- 
ergy spectrum that falls faster than the case without drag 
and £ is the correction to the classical zero drag value of 
the scaling exponent. Note that the major contribution 



to the integrals involved in the theory comes from 
the immediate neighborhood of h = h*, where h\ is the 
value of h at which the minimum of H\{h) occurs. As 
fi — ► 0, h* — ► h and hence £ — ► 0. The scaling expo- 
nent is thus determined by the majority of fluid trajecto- 
ries with stretching rate close to h. On the other hand, 
for jtt ^ 0, h* > h. Therefore the scaling exponent is 
determined by a small number of fluid trajectories that 
experience a stretching rate h\, that is larger than the 
typical stretching rate h, and h\ increases with \x (h\ for 
(j, = 0.1 and fi = 0.2 are 0.47 and 0.57 respectively). The 
reason for this, and in general h* > h as we have shown 
in Section IIII Al is that the presence of drag introduces 
the exponential decaying factor e~^ T in Eq. (|18fl and the 
corresponding integral for E(k) Thus, for a certain 
fixed k (or A), the fluid trajectories that contribute most 
to E(k) (or S2 q (r)) are those with smaller /it, and hence 
larger stretching rate h. 

It should also be clear from the above discussion that, 
in order to get an accurate estimate to £, it is important 
to take into account the fluctuation in the finite-time 
Lyapunov exponent. If we were to neglect the fluctuation 
and take the stretching rate to be h for all trajectories, 
we would have overestimate £ to be 2/j/h. 



2. Vorticity Structure Functions 



The theory predicts that when drag is present, an iner- 
tial range exists in which the vorticity structure functions 
S2q{r) exhibit power-law scaling. The scaling exponent 
C,2 q is given by Eq. (|20[l . In the absence of intermittency 
C,2 q scales linearly with q and (,2qjC,2 — Q, which is plotted 
as the straight dashed line in FIGS.[3» andHJb). For 
/j > 0, Eq. (|20|l predicts that C,2 q varies nonlinearly with 
q, as shown in FIGS. |7Jb) andOJb). This anomalous 
scaling of S2 q (r), which indicates the presence of inter- 
mittency in the system, is verified numerically in both 
cases. 

The anomalous scaling of S2 q (r) is the result of the 
combined effect of drag and non-uniform stretching of 
fluid elements. If all fluid elements have the same stretch- 
ing rate, say h, then e _MT becomes a constant and thus 
&q is proportional to 2q. On the other hand, if /j, = 0, 
regardless of whether the stretching is uniform or not, 
Q 2q = 0. FIGS. Gib) and|Hb) also suggest that S 2q (r) 
shows larger deviation from the simple scaling behavior, 
C29/C2 = 9) as \x increases. 

We remark that the statistical error in the predictions 
of the values of the higher order structure function scaling 
exponents is in general larger. This is because, as already 
mentioned in Section IIII Al h* increases with q, hence 
for large q, (20 depends mostly on a rare number of fluid 
trajectories with large h. 
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FIG. 9: Second order structure functions 52 (r) obtained from 
the energy spectrum E(k) using Eq. 1341 (solid line with cir- 
cles) and from direct numerical simulations (dotted line). 



3. S 2 (r) and E(k) 

We now focus on the second order structure function 
Si (r). With the isotropic assumption, the following re- 
lation between S2M and E(k) can be obtained, 



S 2 (r) 



dk [1 - J {kr)]k 2 E{k) 



(34) 



where Jo is the zeroth order Bessel function of the first 
kind. We compute S%{r) by Eq. J2U) using the E{k) 
shown in FIG. HO The results are plotted as solid lines 
in FIG. |5J the scaling exponents, denoted as (,2,E(k) > are 
then measured and compared with C,2,dns in TABLE [I] 
The S*2 (r) obtained directly from Eq. (|33|) are plotted as 
dashed lines in the same diagram for comparison. Ignor- 
ing viscosity at the small scales and forcing at the large 
scales, we can assume E(k) ~ fc-( 3 +0 for all k. Then it 
follows from Eq. (glj) that, for < £ < 2, S 2 (r) ~ r& for 
all r and C2 = £■ This is consistent with our theory which 
predicts that (2,4/1 = £th when < C2,th < 2. Compar- 
ing (2, dns to £,dns, w e find reasonably good agreements 
but we note that (2, dns is in general slightly larger than 
£dns- We shall see that this small discrepancy is a result 
of the dependence of S2 (r) on the large scales of the flow. 

Due to the effects of forcing and viscosity, E(k) devi- 
ates from the power-law behavior at very small and very 
large k. According to Eq. (|34|) . these deviations will af- 
fect S2(r). This is demonstrated in FIG. ^| in which we 
plot the 5*2 (r) calculated by Eq. (|3"^|) using three different 
E{k). The first one is a pure power law for all k. The sec- 
ond one is the same as the first one except it drops faster 
at large k mimicking the viscous cutoff. The third one 
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FIG. 10: Second order structure functions Sz(r) calculated by 
Eq. ilM II using different E(k): pure power law (cross), power 
law with small scales viscous cutoff (circle) and power law 
with large scales fluctuations (square). The inset shows the 
different E(k) versus k with corresponding labels. 



is constructed from the first one by modifying the first 
five modes to mimic the presence of large scale forcing. 
The range of k used corresponds to a 1024 x 1024 lattice. 
From FIG. it is seen that while the viscous effect is 
negligible except at small r, the large scale forcing can 
have a significant effect on the general shape of S2 (r). 
As a result, S 2 (r) may have a very limited scaling range 
even when E(k) shows a reasonably long one. We be- 
lieve this is generally true for structure functions of any 
order, making it more difficult to accurately measure (29 
than £. We are able to obtain reliable estimates of C29, as 
verified by the agreement between (2, dns and £,dns, by 
employing a 4096 x 4096 lattice and the time-averaging 
process. 



IV. PROBABILITY DISTRIBUTION OF 
VORTICITY DIFFERENCE 

A. Theory 

In this section, we shall derive an expression for the 
probability distribution function P r (dpuj) of the vortic- 
ity difference Spco in terms of the conditional probability 
density function R(t A). In Section llll Al we introduced 
R(t I A) and deduced the scaling of S?ui by doing order of 
magnitude estimations of the integrals in Eq. I|17|). We 
now determine the distribution of Spu) by estimating the 
values of these integrals for different fluid trajectories. 
We shall assume the system is isotropic so that P r {8pw) 
depends on the distance r only, and not on the direction 
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of r . 

Referring to Eq. (|17|l . for t — r < t' < t, we estimate 
that r ■ V/ w [£(*')] ~ {r/L)f u [x{t')] = e~ x f u [x(t')]. For 
< t' < t — t, the separation between the two trajec- 
tories is of order L and the values of f^[x{t')\ along the 
two trajectories become uncorrelated, hence 6f u [x(t')] ~ 
fuj[x(t')]. Approximating f u [x(t')] as varying randomly, 
the values of the first and the second integral in Eq. i|17|) 
are respectively estimated as — h)](e~ x — e _/iT ) 

and (/ a ,//z)e~ M ' r . This suggests that we can treat 5pw as 
a function of two random variables r and Z = f u / fj, as 
follow, 



5fKv(Z, t) j» Z ■ 



(35) 



We then write the probability distribution function of 

6f0j as, 



P r (5?Uj) = P r (5fLu\T)R{T\\)dT (36) 

Jo 

where P r (Spuj\T) is the conditional probability distribu- 
tion function of Spu> given r. From Eq. l|3*K|l and letting 

y ( T ) = (er x + e^ T )/2, we have, 



P r (5pU)\T) 



Pz 



SpU! 

y( T ) 



(37) 



where Pz is the probability density function of Z. Note 
that when r = L, A = r = 0, which implies 8puj\ r= L = Z . 
Therefore, we obtain the result, 



P r (SpLo) 



1 



y( T )' 



-Pl 



y( T ) 



R(t I A) dr 



(38) 



Reference jl3j has previously derived the probability dis- 
tribution function for differences of finite lifetime passive 
scalars advected by a temporally white noise model ve- 
locity field. For the model used in Ref. [lj|] the difference 
probability distribution function at large separation \r\ is 
Gaussian, while this is not the case for our flow. In or- 
der to obtain reasonable agreement between the theory 
Eq. (|38|l and our numerical experiments, it is important 
to account for the non-Gaussian behavior of the proba- 
bility distribution function of the vorticity difference at 
large separation \ f\. In what follows, this will be done by 
using the numerically obtained probability distribution 
function at large \r\ as an input to Eq. l|3*5|) to determine 
the probability distribution function at small |f|. 



B. Comparison of Theory and Numerical Results 

We first compute the probability distribution func- 
tion P r (Xp) of the standardized vorticity difference X?, 
Eq. directly from the numerical solution of Eq. J5J. 
The vorticity field uj(x,t) from t = 41 to t = 75 at every 
1 time unit is used in this computation. The separating 



distance f is taken to be in the ^-direction and is mea- 
sured in the unit of grid size. For fj, = 0.1, the results 
for four different values of r are shown as solid lines in 
FIG. ^2 It is clear that the shape of P r (Xp) changes as 
r varies in the range 4 < r < 512, indicating the system 
is intermittent. P r {Xp) develops exponential tails (e.g., 
r — 64) and then stretched-exponential tails (e.g., r = 4) 
as r decreases. We note that P r (Xp) for all r < 4 col- 
lapse onto each other, and similarly, all the P r (X?) with 
r > 512 have the same shape. Numerical results similar 
to those in FIG. 1111 have also been obtained in Ref. |28j . 
but theory for P r {Xp) was not presented in Ref. psj . 

To apply the theoretical result Eq. l(3*5|) , we first derive 
an expression for R(t | A). To this end, we approximate 
G(h) as a quadratic function of h, 

G(h) « a(h - hf . (39) 

Using the relation Eq. 119f) and the asymptotic form for 
P(h 1 1), Eq. lfTK)l. we obtain 



R(t\X) 



2t £ t 

2 V 7TT V T 



(40) 



To compare the predictions by Eq. (|38|l with numerical 
results, the G(h) for fi — 0.1 is obtained numerically as 
described in Section 111 Rl and fitted to Eq. (|39|) in the 
vicinity of its minimum to obtain the parameters a and 
h. This quadratic fit is a good approximation because 
the integral in Eq. I|38(l is dominated by the region where 
R(t I A) is large, which roughly corresponds to the region 
where P(h 1 1) is maximum. FIG. El shows the R(t \ A) 
for n = 0.1 obtained using the quadratic approximation 
Eq. (EH). We then take P L in Eq. JSHJl to be of the 
form exp(VF) where W is an even sixth degree polyno- 
mial fitted to the numerically computed ln[P r (<5pa;)] for 
r = 1024i and fx = 0.1. With R(t | A) and P L obtained as 
described above, we compute P r (dpuj), and thus P r (Xp), 
for different values of r using Eq. |J5SJ. The results are 
plotted as dashed lines in FIG. ^2 The theoretical pre- 
dictions agree well with the numerical results. We also 
find good agreements between our theory and numerical 
simulations when r is taken to be in the y-direction. Sim- 
ilar results were also obtained for the case fx{k > 6) = 0.2. 



C. Discussion 

According to Eq. I|38|) . P r (dpu>) for a given r can be 
considered as a superposition of many different proba- 
bility distribution functions, each obtained by rescaling 
Pl{S?uj) to a different width using y(r). The amplitude 
of each component is then determined by the coefficient 
R(t\X). Since the statistics between two points sepa- 
rated by distance r ~ L are uncorrelated, the distribu- 
tion Pl(5?uj) is virtually the one-point probability distri- 
bution of the vorticity ui, which is not necessarily Gaus- 
sian and is closely related to the statistics of the source 
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FIG. 11: For /j, = 0.1, the probability distribution function P r (Xp) of the standardized vorticity difference X? obtained from 
direct numerical simulation (solid lines) and from Eq. (1381 (dashed lines) . The separating distance r is in the ^-direction and 
measured in the unit of grid size. 



as implied by the definition of Z. Hence, Eq. I|38|) 
relates the distribution of the vorticity difference to the 
one-point statistics of the source. 

As can be seen in FIG. El R{ T I A) has a very sharp 
peak when A is small (r is large). Thus, only a small range 
of values of r contributes to the integral Eq. I|38|) . This 
gives the expected result that P r {8poj) has similar shape 
as Pl{8fu) when r is large. As A increases (r decreases), 
the width of R(t | A) increases. According to Eq. l|55)l. 
P r (8puj) is now composed of a large number of rescaled 
Pl{5poj) with a broad range of width. This gives rise to 
the broader-than- Gaussian tails in P r {8puj) for small r. 

Like the anomalous scaling of the vorticity structure 
functions, the dependence of the shape of P r (8puj) on r 
is a result of the presence of drag and the non-uniform 
stretching of fluid elements. If h = h for all fluid tra- 
jectories, then R(t | A) = S(t — X/h), and from Eq. ffiBfr. 
P r (8?u)) has the same form for all r, indicating a self- 
similar flow. On the other hand, if ji = 0, y(r) becomes 
independent of r, which also implies P r (8?Lo) is indepen- 
dent of r. In both situations, P r (8pLu) will have the same 
shape as Pl(Spuj). 



We remark that while the width of the distribution 
P r (8puj) in Eq. (f^H|) depends precisely on the choice of 
the value for the multiplicative constant on the right hand 
side of Eq. (|35|l . our results for P r (X?) are independent 
of this value. 



V. MULTIFRACTAL FORMULATION 
A. Theory 

The local viscous energy dissipation rate per unit mass 
e and its global average value (e) play important roles in 
the phenomenology of three-dimensional turbulence |30| . 
It is now well known that e shows intermittent spatial 
fluctuations which can be described by the multifractal 
formulation |31| . Using the measure p'(r) = £ r /£, where 
£ r is the total dissipation in a volume of linear dimension 
r and £ is the total dissipation in the whole domain, 
the Renyi dimension spectrum D q and the singularity 
spectrum f(a) have been measured experimentally [321 
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FIG. 12: For /i = 0.1, the conditional probability density 
function R(t | A) given by Eq. (a = 1.55, h = 0.26) for 
different values of A. 



Intermittency in three-dimensional turbulence also 
manifests itself as anomalous scaling in the velocity struc- 
ture functions S' 3q (r) defined as 



S' 3q (r) = ((5,v)^)^r^ 



(41) 



where v is the component of the velocity vector in the di- 
rection of rand 5?v = v(x+r)—v(x). From Kolmogorov's 
hypotheses in his 1941 paper (3^|, which ignore the inter- 
mittency of e, one arrives at the result ^ q = q. Exper- 
iments show that £' 3q is a nonlinear function of q. This 
anomalous scaling of S' 3q (r) is believed to be related to 
the intermittency of e. Kolmogorov's 'refined similarity 
hypothesis' in his 1962 paper |35| gives the connection 
between intermittency in velocity difference and inter- 
mittency in e. The refined similarity hypothesis states 
that at very high Reynolds numbers, there is an inertial 
range of r in which the conditional moments of Spv scales 
as follows 



((5fv)' 1 | e r ) ~ (re r ) 



q/3 



(42) 



where e r is the average of e over a volume of linear dimen- 
sion r. Eq. 142(1 implies {(5pv) q ) ~ {(re r ) q ^ 3 ) which gives 
the following relation between the D q based on p'(r) and 



D„ 



a q - <zC 3 



(43) 



Kolmogorov's fourth-fifths la w g ives £3 = 1 exactly |31j . 
Eq. H43|) can also be derived [33] from the relation 



(M 3 



(44) 



We note that Eq. follows from Eq. |@2J|. As shown 
below, there is a relation analogous to Eq. (|43J) for the en- 
strophy cascade of two-dimensional turbulence with drag. 

For the enstrophy cascade regime in two-dimensional 
turbulence, a central quantity to the phenomenology Q 
in this regime is the viscous enstrophy dissipation 77 given 
by Eq. Q, and the relevant measure is Pi(e) defined in 
Eq. @ . We have already seen in Section IIIII that in 
the presence of drag, the vorticity structure functions 
scale anomalously with scaling exponents £2,3 given by 
Eq. (|20l) . We now derive a relation between <^ 2q and the 
D q based on Pi(e). 

Consider the following quantity 



(45) 



which by the definition of D qi Eq. (|10(1 . scales like 

h{q,e)~e^ D * (46) 



for some range of e. Assume there exists a scaling range 
extending from the system scale L ~ kj 1 down to the 

dissipative scale r d ~ k^ 1 such that both the scaling 
relations Eq. (@J and Eq. hold. At the dissipative 
scale, due to the action of viscosity, the vorticity field tu 
becomes smooth, thus we have the following relations, 



K{r d ) 



\Vu\ 2 dx ~ r 2 \Vuj(x)\ 2 , f£%) (47) 



IVwl 



\r\=r d . (48) 



Then by putting e = in Eq. © and let |r| = and 

xi e Hi(rd), we get 



Pi(r d ) 



\SpL0(Xi)\ 2 



( r d) 2 J n \SrUl\ 2 dx 



(49) 



Substituting Eq. (|4*9"|l in Eq. (|4*5|l . we obtain the scaling 
of h(q, e), 



Y J \5Mx l )\ 2q 



h(q,r d ) 



(r d ) 2 i 



\5fU)\ 2 dx 



-) (l^-l 29 ) 
r d 



29 



(\Spuj 



(r d ) 



2<?-2+C 2 ,-9C2 



(50) 



Comparing Eq. l(4^|) to Eq. JSUJ), we get the principal 
result of this section. 



D q = 2 



(2q ~ q& 

9-1 



(51) 
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which can be regarded as analogous to Eq. (|43() . 

We mention that Eq. II. 'ill can be derived in an anal- 
ogous manner using Eq. (14411 . From Eq. (|51|l . in two- 
dimensional turbulence, the measure pi is multilractal 
when the vorticity structure lunctions exhibit anomalous 
scaling. Hence, in the presence of drag, we expect the 
measure based on the squared vorticity gradient |Vcj 
to show multifractal structures. 



2 





B. Comparison of Theory and Numerical Results 



The multifractal structure of |Vw| 2 is most readily 
visualized in snapshots of |Vo;| 2 from our simulations. 
Since |Vw| 2 grows at widely varying exponential rates, 
only a few points would be visible if |Vw| 2 was plotted 
directly using a linear scale. Therefore, we plot the fol- 
lowing quantity instead |27j . 



M(x) 



(52) 



where the set Kg contains those lattice points Xi for which 
|Vw(afi)| 2 < |Vu;(af)| 2 , and we sum over all lattice points 
in the denominator. By definition, < Ai(x) < 1. 
FIG. El shows the results for /i = 0.1 and // = 0.2. 
Filament structures can clearly be seen for both cases, 
showing that the measure Pi concentrates in a very small 
area. This is particularly clear in the case fj, = 0.2. 

To quantify the multifractal nature of Pi, we now cal- 
culate its Renyi dimension spectrum D q . We employ the 
box-counting method to estimate D q . Using box size e/L 
ranging from 2~ 12 to 2 _1 , we compute the instantaneous 
Ii(q, e) from t = 41 to t = 75 at every 1 time unit using 
the numerical solution of Eq. For q ^ 1, we then 
make log-log plot of the time-average of [Ii/(q — 1)] ver- 
sus e/L. These are shown in FIGS. Et a ) and llftf a). For 
q = 1, we take the limit q — > 1 in Eq. (|46|1 to obtain 



He) 



where 



He) 



Pi(e) log 2 [K(e)] 



(53) 



and for q = 1 in FIGS. Hlf a) and I15f a). the time aver- 
age of /2(e) is plotted against log 2 (e/L). According to 
Eq. (|46() , these plots will show a linear region with slope 
equals D q . All curves in FIGS. I14f a) and I15f a) show 
slightly undulating behavior which introduces uncertain- 
ties in the determination of D q . The estimated D q at 
different values of q are shown as circles with error bars 
in FIGS, n^r b) and llSf b). The error bars correspond to 
the variability of the D q observed at different moments 
in time. The dotted lines in the figures are fourth degree 
polynomials fitted to the circles. We also compute D q 
using Eq. (|5T)l . To this end, we fit the curves of C,i q jC,i 
versus q in FIG. □» and FIG. Efb) with the following 
polynomial, where d — 3 for the circles and d = 5 for the 




FIG. 13: Snapshots of the scaled squared vorticity gradient 
|Voj| 2 at t — 61 for the case /i = 0.1 (upper) and at t — 65 
for the case /1 = 0.2 (lower). Light areas are regions of large 
values, and dark areas are regions of small values. 



C2 



(54) 



By Eq. J5TJ, D q is then given by 



D q = 2 + ( 2 J2anq(q-l) n - 1 

n=l 



(55) 
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FIG. 14: For the case of fi — 0.1: (a) Ii(q, e) for q between 
0.1 and 2.0 (7 2 (e) is plotted for q = 1.0). The dotted lines 
are linear fits in the scaling region, (b) D q computed using 
numerical solution of Eq. ||5J (circle with error bar) and its 
fourth degree polynomial fit (dotted line). D q predicted by 
the theory Eq. 15111 when C,2 q obtained from numerical simula- 
tions are used (square) and when <^2 q calculated from Eq. I20H 
are used (diamond). 
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FIG. 15: For the case of fi — 0.2: (a) Ii(q,e) for q between 
0.1 and 2.0 (I 2 (e) is plotted for q = 1.0). The dotted lines 
are linear fits in the scaling region, (b) D q computed using 
numerical solution of Eq. J5J (circle with error bar) and its 
fourth degree polynomial fit (dotted line). D q predicted by 
the theory Eq. 15111 when £29 obtained from numerical simula- 
tions are used (square) and when C : 2 q calculated from Eq. I20H 
are used (diamond). 



In FIGS. [TUb) andEtb), we plot Eq. (jSEJ using C, 2q 
obtained from numerical simulations, as well as £2? cal- 
culated from our theory. The results are shown as solid 
lines labeled with squares and diamonds, respectively. 
Despite the fact that there are discrepancies between the 
D q obtained by the various methods, they all show the 
same trend and clearly indicates that pi is multifractal 



(i.e., D q varies with q). 

We now generate the singularity spectrum f(a) 
by Legendre transforming the D q curves shown in 
FIGS. fHTT b) and !15f b). In particular, for each value of 
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a 

FIG. 16: For the case of fi = 0.1: The /(a) curves generated 
by Legendre transforming the D q curves in FIG. mf bl with 
the corresponding labels (circle, square and diamond), /(a) 
at different values of a obtained by the canonical method is 
shown as well (cross). 

q, we have 

a(q) = ^-[{q-l)D q ] , (56) 
f[a(q)} = aq-(q-l)D q . (57) 

In the absence of intermittency D q is independent of q 
and f(a) is defined only at a — D q . Thus the consis- 
tent determination of f(ct) over some range of a indicates 
intermittency. FIG. [ED and FIG. IT7I show the f(a) ob- 
tained for jj, — 0.1 and /i = 0.2 respectively using the D q 
obtained by the three methods in FIG. I14f b1 and !15f b) 
(circle, square and diamond symbols). The results for 
f{a) are seen to agree very well with each other in spite 
of the difference between the D q determinations. Thus, 
we find that f(a) gives a more consistent measure of 
intermittency across different methods of determination 
than does D q . We believe that the disagreement seen 
in FIGS. H4T b) and !15f b) between the different methods 
for determining D q is not significant in view the limited 
amount of scaling range available. 

We can also determine f(a) directly from the numer- 
ical solution of Eq. (J2J). Following Ref. [32ll. we use the 
canonical method developed in Ref. j3|j to determine 
/(a). Accordingly, we construct the normalized q th order 
measures rrii{q,e) in box IZi(e) as follow, 

'"•"• f, = Efp' (58) 

Then, the mean singularity index a and the correspond- 
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FIG. 17: For the case of /i = 0.2: The f(a) curves generated 
by Legendre transforming the D q curves in FIG. I15f b'l with 
the corresponding labels (circle, square and diamond). f(a) 
at different values of a obtained by the canonical method is 
shown as well (cross). 

ing f(a) is given by 

a(q) = rrnoLi = lim lim , % — %J^1. (59) 
t-r 1 e^o^o loge 

f\a(q)] = lim lim ^™* lo Z m * . (60) 
J L WJ e^^o loge v ' 

One of the advantages of the canonical method is the 
absence of finite-size effect due to logarithmic prefac- 
tors . To estimate a, we plot the time average of 
the quantity J2i m i l°g2 Pi versus log 2 (e/L) and measure 
the slopes in the scaling region. Similarly, when we plot 
the time average of ^ i TOilog 2 rrii versus log 2 (e/L), the 
slopes in the scaling range give the values of /(a). For 
the case q = 1, a(l) = /[a(l)] is obtained by measuring 
the slope of the curve of the time-averaged I^{l f e) versus 
log 2 (e/L). The values of /(a) at different values of a 
obtained by this scheme are shown as crosses in FIG. 1161 
and FIG. El for fj, = 0.1 and /i = 0.2 respectively. They 
are in excellent agreement with the f{a) generated by 
Legendre transforming the D q curve obtained from the 
numerical solution of Eq. © . 

C. Discussion 

Kolmogorov introduced the refined similarity hypoth- 
esis (RSH) Eq. 142(1 to take into account the spatial fluc- 
tuations of e in three-dimensional turbulence. The re- 
lation Eq. (|43|l is a direct consequence of the RSH. In 
two-dimensional turbulence, the relevant quantity is the 
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local rate of viscous enstrophy dissipation rj. We have 
already seen that the measure pi is multifractal in the 
presence of drag, indicating the intermittent nature of rj. 
Following Kolmogorov's ideas, we consider the average of 
77 over an area lZ(r) of linear dimension r: 



v 

^ = -2 



Vwl 2 dx . 



(61) 



-R{r) 



Analogous to the Kolmogorov's RSH, we propose that at 
high Reynolds number, there is an inertial range of r in 
which 



(\6M 2q \r) r ) = C 2q (r^r, r y 



(62) 



where C 2q are constants. Eq. (|62[1 implies (\8fuj\ 2q ) 
{(Vr) q )r qC2 and hence {(r]r) q ) ~ r 7 ' with 



Iq = (2q ~ qC.2 ■ 



(63) 



An expression analogous to Eq. (|63[1 has been proven for 
Kraichnan's model of passive scalar advection [3?], l38| . 
We note that Eq. follows if 



Vr 



(64) 



which is the two-dimensional counterpart of Eq. I|44|) . It 
is straightforward to show that the hypothesis Eq. 162(1 
implies Eq. 1)510. Denoting 77,. in box lZi(r) by rjr l \ we 
have 



Pi(r) 



2 (i) 

r r]r 



(65) 



Since X^I 7 ? 



(Oi 



r 2 {{ r lr) q ), we get I(q,r) ~ r 



2q-2+7, 



from which Eq. I|51|) immediately follows. 

Finally, we remark that if the measure pi is defined 
with a parameter n as follow p7|. 



S n \Vu\ n dg 



(66) 



then the corresponding formula for the D q based on this 
measure is 



D q = 2 



Cnq qCn 



q 



1 



(67) 



VI. CONCLUSION 

We have studied the enstrophy cascade regime of two- 
dimensional turbulence with linear drag. A previous the- 
ory for the power law exponent of the energy wavenum- 
bcr spectrum is verified by direct numerical computa- 
tion using a 4096 x 4096 lattice. We also calculate the 
vorticity structure functions numerically and show that 
they exhibit anomalous scaling in the presence of drag. 
The values of the structure function scaling exponents £29 
are measured and found to agree with the prediction by 
previous theory. We then compute the probability dis- 
tribution function P r (X?) of the standardized vorticity 
difference X? and find that P r (Xp) develops exponential 
and stretched-exponential tails at small values of r. The 
theoretical expression for P r (X?) is shown to give pre- 
dictions that agree well with the numerical results for a 
wide range of r. A measure based on the local viscous en- 
strophy dissipation rate 77 is studied in terms of its Renyi 
dimension spectrum D q and singularity spectrum f(a), 
and is found to be multifractal. The intermittency in 77 
is connected to the intermittency in vorticity difference 
by a two-dimensional analog of the refined similarity hy- 
pothesis, and we derive a formula that relates D q to ^ 2q . 
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